(*prefix=(SetDirectory[NotebookDirectory[]]//ParentDirectory//ParentDirectory//ParentDirectory//ParentDirectory//ParentDirectory)<>"/Articles/NeutrinoCooledStar/calc/"*)

prefix="https://th.if.uj.edu.pl/~odrzywolek/psns/data/";

Dla masywnych gwiazd na ostatnich etapach ewolucji, lata przed wybuchem supernowej, diagram HR staje się bezużyteczny. Wszystko co istotne dzieje się w jądrze gwiazdy, a produkowane promieniowanie ma formę neutrin. W takiej sytuacji diagram prezentujący jasność absolutną i temperaturę powierzchni należy zastąpić  jasnością neutrinową i średnią energią neutrin. Sytuacja staje się nieco bardziej skomplikowana, gdyż mamy efektywnie  trzy rodzaje neutrin:  (neutrina elektronowe)  (antyneutrina elektronowe) i zbiorczo 4 pozostałe zapachy neutrin (, , , ).  Co więcej, neutrina podlegają oscylacjom (w materii), co może powodować zmiany ich zapachu w trakcie podróży ze źródła do detektora.

OMKs15tot1=Import[prefix<>"s15_1.tot.dat","Table"];
OMKs15tot2=Import[prefix<>"s15_2.tot.dat","Table"];
OMKs15tot3=Import[prefix<>"s15_3.tot.dat","Table"];
nums15=Dimensions[OMKs15tot1][[1]]

Lν1s15=Table[{1/OMKs15tot1[[ii]][[1]],OMKs15tot1[[ii]][[2]]},{ii,nums15}];
Lν2s15=Table[{1/OMKs15tot2[[ii]][[1]],OMKs15tot2[[ii]][[2]]},{ii,nums15}];
Lν3s15=Table[{1/OMKs15tot3[[ii]][[1]],OMKs15tot3[[ii]][[2]]},{ii,nums15}];

Nν1s15=Table[{OMKs15tot1[[ii]][[1]],OMKs15tot1[[ii]][[3]]},{ii,1,nums15}];
Nν2s15=Table[{OMKs15tot2[[ii]][[1]],OMKs15tot2[[ii]][[3]]},{ii,1,nums15}];
Nν3s15=Table[{OMKs15tot3[[ii]][[1]],OMKs15tot3[[ii]][[3]]},{ii,1,nums15}];

timeticksINV={{1/(36525*86400),"100 yr"},{1/(3652.5*86400),"10 yr"},{1/(365.25*86400),"1 yr"},{1/(100*86400),"100d"},{1/(10*86400),"10d"},{1/86400,"1d"},{1/(86400/8),"3h"},{1/(86400/24),"1h"},{1/(86400/24/6),"10 min"}};

timeticks={{(36525*86400),"100 yr"},{(3652.5*86400),"10 yr"},{(365.25*86400),"1 yr"},{(100*86400),"100d"},{(10*86400),"10d"},{86400,"1d"},{(86400/8),"3h"},{(86400/24),"1h"},{(86400/24/6),"10 min"}};

SetOptions[ListLogLogPlot,
Joined->True,
Frame->True,
AspectRatio->1/2,
FrameTicks->{{ Automatic,Automatic},{timeticks,None}},FrameLabel->{{"log \!\(\*SubscriptBox[\(L\), \(ν\)]\) [\!\(erg\*SuperscriptBox[\(s\), \(-1\)]\)]",None},{"Time B.C.",None}},
BaseStyle->{FontSize->14},PlotStyle->{{Blue,Dashed,Thickness[0.004]},{Red,Thickness[0.004]},{Black,Dotted,Thickness[0.004]},{Black,Dotted,Thickness[0.004]}},Joined->True,
PlotRange->{{600,36600*86400},{10^46,10^54}},
ImageSize->600,
FrameStyle->Thickness[0.003]];

s15Flux=ListLogLogPlot[{Nν1s15,Nν2s15,Nν3s15},PlotLabel->"15 (Subscript[M, Sun]) (s15)",PlotLegends->{"(Subscript[ν, e])","(Subscript[(OverBar[ν]), e])","(Subscript[ν, μ]),(Subscript[(OverBar[ν]), μ])"}]

Eν1s15=Table[{OMKs15tot1[[ii]][[1]],OMKs15tot1[[ii]][[4]]},{ii,1,nums15}];
Eν2s15=Table[{OMKs15tot2[[ii]][[1]],OMKs15tot2[[ii]][[4]]},{ii,1,nums15}];
Eν3s15=Table[{OMKs15tot3[[ii]][[1]],OMKs15tot3[[ii]][[4]]},{ii,1,nums15}];

kpc=3.08568025*10^21

yminorticks2=Flatten[Table[Table[{N[4 Pi kpc^2 xx 10^yy],""},{xx,2,9,2}],{yy,1,10}] ,1];
ymajorticks2=Table[{4 Pi kpc^2 10^xx,ToString[xx],{0.02,0}},{xx,1,10}];
yticks2=Union[ymajorticks2,yminorticks2];

SetOptions[ListLogLogPlot,
Joined->True,
Frame->True,
AspectRatio->1/2,
FrameTicks->{{ Automatic,Automatic},{timeticksINV,None}},FrameLabel->{{"log \!\(\*SubscriptBox[\(L\), \(ν\)]\) [\!\(erg\*SuperscriptBox[\(s\), \(-1\)]\)]",None},{"Time B.C.",None}},
BaseStyle->{FontSize->14},PlotStyle->{{Blue,Dashed,Thickness[0.004]},{Red,Thickness[0.004]},{Black,Dotted,Thickness[0.004]},{Black,Dotted,Thickness[0.004]}},Joined->True,
PlotRange->{{10^-9,10^-1},{10^40,10^50}},
ImageSize->600,
FrameStyle->Thickness[0.003]];

s15noOSC=ListLogLogPlot[{Lν1s15,Lν2s15,Lν3s15},PlotLegends->{"(Subscript[ν, e])","(Subscript[(OverBar[ν]), e])","(Subscript[ν, μ]),(Subscript[(OverBar[ν]), μ])"},PlotLabel->"No oscillations (s15 model)"]

Wprowadzamy poprawkę na oscylacje neutrin w materii.

p1NORM=0.0234

p2NORM=0.676

timeBCs15=Lν1s15[[;;,1]];

Lν1oscNORM={timeBCs15,p1NORM Lν1s15[[;;,2]] +(1-p1NORM)Lν3s15[[;;,2]]}//Transpose;

Lν2oscNORM = {timeBCs15, p2NORM Lν2s15[[;;,2]] +(1-p2NORM)Lν3s15[[;;,2]]}//Transpose;

Lν3oscNORM={timeBCs15, (1-p1NORM) Lν1s15[[;;,2]] +p1NORM Lν3s15[[;;,2]]}//Transpose;

Lν4oscNORM={timeBCs15, (1-p2NORM) Lν1s15[[;;,2]] +p2NORM Lν3s15[[;;,2]]}//Transpose;

s15NORM=ListLogLogPlot[{Lν1oscNORM,Lν2oscNORM,Lν3oscNORM,Lν4oscNORM},PlotLegends->{"(Subscript[ν, e])","(Subscript[(OverBar[ν]), e])","(Subscript[ν, μ])","(Subscript[(OverBar[ν]), μ])"},PlotLabel->"Normal ν mass hierarchy (s15 model)"]

p1INV=0.3

p2INV=0.0234

Lν1oscINV={timeBCs15,p1INV Lν1s15[[;;,2]] +(1-p1INV)Lν3s15[[;;,2]]}//Transpose;

Lν2oscINV = {timeBCs15, p2INV Lν2s15[[;;,2]] +(1-p1INV)Lν3s15[[;;,2]]}//Transpose;

Lν3oscINV={timeBCs15, (1-p1INV) Lν1s15[[;;,2]] +p1INV Lν3s15[[;;,2]]}//Transpose;

Lν4oscINV={timeBCs15, (1-p2INV) Lν1s15[[;;,2]] +p2INV Lν3s15[[;;,2]]}//Transpose;

s15INV=ListLogLogPlot[{Lν1oscINV,Lν2oscINV,Lν3oscINV,Lν4oscINV},PlotLegends->{"(Subscript[ν, e])","(Subscript[(OverBar[ν]), e])","(Subscript[ν, μ])","(Subscript[(OverBar[ν]), μ])"},PlotLabel->"Inverted ν mass hierarchy (s15 model)"]

s15nue=ListLogLogPlot[{Lν1s15, Lν1oscINV,Lν1oscNORM},PlotLegends->{"No oscillations", "Inverted", "Normal"},PlotLabel->"Electron neutrino luminosity (s15 model)",PlotStyle->{{Black,Thickness[0.004]},{Black,DotDashed, Thickness[0.004]},{Black,Dashed,Thickness[0.004]}}]

s15nuebar=ListLogLogPlot[{Lν2s15, Lν2oscINV,Lν2oscNORM},PlotLegends->{"No oscillations", "Inverted", "Normal"},PlotLabel->"Electron anti-neutrino luminosity (s15 model)",PlotStyle->{{Black,Thickness[0.004]},{Black,DotDashed, Thickness[0.004]},{Black,Dashed,Thickness[0.004]}}]

s15shading=LogLinearPlot[5(UnitStep[(t-25*60)(1.7*3600-t)]+UnitStep[(t-12.9*3600)(8.3*86400-t)]+UnitStep[(t-16.75*86400)(2.6*365.25*86400-t)]),{t,600,100*365*86400},PlotRange->{{600,36600*86400},{0,5}},Frame->True,Filling->Bottom,FillingStyle->Opacity[0.2]]

s15shading2=LogLogPlot[10^55(UnitStep[(t-25*60)(1.7*3600-t)]+UnitStep[(t-12.9*3600)(8.3*86400-t)]+UnitStep[(t-16.75*86400)(2.6*365.25*86400-t)]),{t,600,100*365*86400},PlotRange->{{600,36600*86400},{10^46,10^54}},Frame->True,Filling->Bottom,FillingStyle->Opacity[0.2]]

SetOptions[ListLogLinearPlot,Joined->True,Frame->True,AspectRatio->(9/16),FrameTicks->{Automatic,Automatic, timeticks,None},FrameLabel->{"Time B.C. [seconds]","Avg. ν energy [MeV]",None,None},BaseStyle->{FontFamily->"Arial",FontSize->14},PlotRange->{{600,36600*86400},{0,5}},PlotStyle->{{Black,Dashed},Black,{Black,Dotted}},ImageSize->450];

s15Enu=ListLogLinearPlot[{Eν1s15,Eν2s15,Eν3s15},PlotLabel->"15 (Subscript[M, Sun]) (s15)"]

Show[s15Enu,s15shading,Epilog->{Inset[Rotate["Shell Si\nburning",1],{Log[3000],0.5},BaseStyle->{FontFamily->"Arial",FontSize->10}],Inset[Rotate["Core Si\nburning",1],{Log[2*10^5],0.5},BaseStyle->{FontFamily->"Arial",FontSize->10}],Inset[Rotate["Core/shell Oxygen\nburning",1],{Log[10^7],2.4},BaseStyle->{FontFamily->"Arial",FontSize->10}]},ImageSize->450]

Show[s15Flux,s15shading2,ImageSize->450,Epilog->{Inset["Shell Si\nburning",{Log[3000],Log[10^47]},BaseStyle->{FontFamily->"Arial",FontSize->10}],Inset["Core Si\nburning",{Log[2*10^5],Log[10^47]},BaseStyle->{FontFamily->"Arial",FontSize->10}],Inset["Core/shell Oxygen\nburning",{Log[10^7],Log[10^47]},BaseStyle->{FontFamily->"Arial",FontSize->10}]},FrameStyle->Thickness[0.002]]

Ponieważ na ostatnich etapach spalania jądrowego (C, O, Si burning) pre-supernowych procesy w jądrze odbywają się zbyt szybko aby uwidocznić się na klasycznym diagramie HR, potrzebujemy jego wersji neutrinowej. Na osi poziomej pozostaje jasność, ale tym razem neutrinowa, która jest o kilka-kilkanaście rzędów wielkości większa od fotonowej na tym etapie ewolucji. Zamiast temperatury (typu widmowego) na osi poziomej umieszczamy średnią energię neutrina. Jest ona zależna od temperatury, w sposób bardzo podobny jak dla gazu fotonowego (ciała doskonale czarnego). Dla neutrin  istotna jest jednak degeneracja materii (neutrina, elektrony i nukleony są fermionami), co powoduje, że średnią energię podnosi potencjał chemiczny μ. W  odróżnieniu od fotonów, dla których oczywiście μ=0. Ponadto, energia neutrina ściśle powiązana jest z procesami jego produkcji, co łatwo pozwala zidentyfikować na diagramie różne fazy ewolucji. Np. gdy średnia energia neutrina przekracza energię spoczynkową elektronu, 0.511 MeV/c^2, oznacza to pojawienie się procesu anihilacji par (SuperPlus[e])+(SuperMinus[e]) na pary neutrino-antyneutrino. Co więcej, procesy detekcji neutrin są często odwrotne do tych je produkujących, i diagram OMK pozwala na odczytanie możliwości ich detekcji.

nuHR1s15=Table[{OMKs15tot1[[ii]][[4]],OMKs15tot1[[ii]][[2]]},{ii,1,nums15}];
nuHR2s15=Table[{OMKs15tot2[[ii]][[4]],OMKs15tot2[[ii]][[2]]},{ii,1,nums15}];
nuHR3s15=Table[{OMKs15tot3[[ii]][[4]],OMKs15tot3[[ii]][[2]]},{ii,1,nums15}];

Diagram OMK. Są trzy linie, bo mamy 3 istotnie różne neutrina: (Subscript[ν, e]) (neutrino elektronowe)  (antyneutrino elektronowe)  (neutrino “mionowe” - w astrofizyce pozostałe 4 zapachy neutrin , , ,  ,,binujemy’’ do jednego, a neutrino mionowe jest ich reprezentantem). Neutrina mionowe dają ,,czysty’’ sygnał procesów termicznych: anihilacji, rozpadu plazmonu, fotoprodukcji czy brebstrahhlungu. Neutrina elektronowe, szczególnie zaczynając od spalania Si, są produkowane w procesach jądrowych (rozpady β), głównie wychwycie elektronu. Proces ten zachodzi także, ale znacznie słabiej (z powodu braku pozytonów), dla antyneutrin elektronowych.

ListLogPlot[{Legended[nuHR3s15,"(Subscript[ν, μ])"],Legended[nuHR2s15,"(OverBar[(Subscript[ν, e])])"],Legended[nuHR1s15,"(Subscript[ν, e])"]},Joined->False,InterpolationOrder->1,PlotRange->{{0.0,4},{10^30,10^50}},Frame->True,AspectRatio->4/3,GridLines->{{0.511},None}]